Quantum control using sequences of simple control pulses 
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Structured decompositions of a desired unitary evolution operator are employed to derive control 
schemes that achieve certain control objectives for finite-level quantum systems using only sequences 
of simple control pulses such as square-waves with finite rise and decay times or Gaussian wavepack- 
ets. The technique is applied to find control schemes that achieve population transfers for pure-state 
systems, complete inversions of the ensemble populations for mixed-state systems, create arbitrary 
superposition states and optimize the ensemble average of observables. 
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I. INTRODUCTION 

The ability to control quantum systems is an essen- 
tial prerequisite for many present and potential future 
applications involving the manipulation of atomic and 
molecular quantum states 0. Some of the important 
applications are quantum state engineering Q , control of 
chemical reactions || and laser cooling of internal molec- 
ular degrees of freedom ||, pj. Quantum control theory 
may also reveal new ways to solve problems crucial to 
quantum computing such as initializing quantum regis- 
ters H or building robust quantum memory 0. 

Due to the wide range of applications, the immediate 
aims of quantum control may vary. However, the control 
objective can usually be classified as one of the following: 

1 . to steer the system from its initial state to a target 
state with desired properties, 

2. to maximize (or minimize) the expectation value 
(ensemble average) of a selected observable, or 

3. to achieve a certain evolution of the system. 

Despite their apparent dissimilarity, these control objec- 
tives are closely related. Indeed, (1) can be considered 
a special case of (2) in which the observable is the pro- 
jector onto the subspace spanned by the target state and 
(2) is a special case of (3) where we desire to find an evo- 
lution operator that maximizes the expectation value of 
the selected observable either at a specific target time or 



at some time in the future. Thus, one of the central prob- 
lems of quantum control is to achieve a desired evolution 
of the system by applying external control fields and the 
primary challenge is to find control pulses (or sequences 
of such pulses) that are feasible from a practical point of 
view and effective in realizing the control objective. 

Many of the control strategies for quantum systems 
that have been proposed to date rely on numerical 
methods to solve the control or optimization problem 
|| @ H H Q. In this paper, we explore 
an alternative approach to quantum control, which is 
based on explicit generation of unitary operators using 
Lie group decompositions The technique uses the 

rotating wave approximation (RWA) and employs fre- 
quency discrimination or other atomic selection rules to 
address transitions individually. Although these assump- 
tions seem to preclude the application of the technique 
to some model systems of interest, e.g., systems with 
equally spaced or almost equally spaced energy levels, it 
should nevertheless be applicable to many quantum sys- 
tems. A significant advantage of the approach we pursue 
in this paper is the flexibility of the control pulse shapes, 
which makes it possible to choose control pulses that are 
reasonably simple to create in the laboratory, such as 
square-wave pulses with finite rise and decay times or 
Gaussian wavepackets. 



II. MATHEMATICAL FORMULATION OF THE 
CONTROL PROBLEM 



'Electronic address 
t Electronic address 
* Electronic address 
§ Electronic address 



*.p\schirmer@open.aciik 




n. ac.uk 


vislrOmtdallas.edu 




irabitz@chcmvax.princeton.edr 



We consider a quantum system with a discrete energy 
spectrum whose energy levels E n are assumed to be finite, 
non-degenerate and ordered in an increasing sequence, 



E x < E 2 < ■ ■■ < E_ 



N, 



(1) 



2 



E =2.275 

4 



E 3 =1.875 



E 2 =1.275 



E 1 =0.475 



H 3 =0.4, d 3 =3 



1/2 



H 2 =0.6, d 2 =2 



1/2 



I 



1^=0.8, d 1= 1 



FIG. 1: Energy-level and transition diagram for a four-level 
Morse oscillator with anharmonicity a = 0.1. 



where N < oo is the dimension of the Hilbert space of 
pure states of the system. Note that non-degeneracy and 
the particular ordering of the energy levels are assumed 
primarily for simplicity. These requirements can often be 
relaxed if selection rules are available to distinguish de- 
generate energy levels and different transitions that have 
the same transition frequency. Moreover, non-degeneracy 
of the energy levels can often be assured by considering 
suitable subspaces of the Hilbert space of pure states. 
Our model is strictly valid only for dissipation-free sin- 
gle atoms or optically thin ensembles of homogeneously 
broadened atoms. However, by choosing appropriately 
short control times (with respect to the lifetimes of the 
excited states) jl6| , our idealized model should be appli- 
cable to a variety of quantum systems including many 
atomic and molecular iV-state systems, Rydberg atoms 
or particles in (anharmonic) quantum wells. 

Although the results presented in this paper do not 
depend on a specific model, for clarity, we shall illustrate 
our results using a four-level Morse oscillator with energy 
levels 
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for 1 < n < 4 and transition dipole moments 



P12 



\/n, 1 < n < 3, 



(2) 



(3) 



where ujq andpi2 are constants representing the oscillator 
frequency and the transition dipole moment of the transi- 
tion |1) —>■ 1 2), respectively. The parameter a determines 
the anharmonicity of the system, which we arbitrarily set 
to 0.1. The energy-level and transition diagram for this 
system is shown in figure [l| 

The free evolution of the system is governed by the 



internal Hamiltonian, whose spectral representation is 
//„ >]E n \n)(n\, (4) 
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where {\n} : n = 1, . . . , N} is a complete set of orthonor- 
mal eigenstates that satisfy the stationary Schrodinger 
equation 



H \n) =E n \n), l<n<N. 



(5) 



The application of external control fields perturbs the 
system and gives rise to a new Hamiltonian H = Hq+Hj, 
where Hi is an interaction term. For a control field of 
the form 



f m {t) = 2A m (t) cos(/z m t + (j> m ) 



A m (i) 



(6) 



which is resonant with the transition frequency fj, m = 
E m +i —E m and drives only the transition \m) — ► |m + l), 
using the rotating wave approximation (RWA) leads to 
the following interaction term 

H m (frn) = A m (t)e^ > d m \ Ul) (m + 1 1 

+A m (t)e-^ t+ ^d m \m+l)(m\, (7) 

where d m is the transition dipole moment for the transi- 
tion | to) — > |to + 1), 4> m is the initial phase and 2A m {t) 
is the envelope of the field, which must be slowly varying 
with respect to /x" 1 . It should be noted that the validity 
of this approximation for the interaction with the field f m 
depends on one's ability to individually address the mth 
transition, either by use of selection rules or frequency 
discrimination. 

For frequency-selective control the transition frequen- 
cies n m must be sufficiently distinct to assure that each 
field is resonant only with a single transition and that off- 
resonant effects on other transitions are small. To satisfy 
this requirement, we note that if a monochromatic field 
of constant amplitude 2A m and frequency fj, m is applied 
then off-resonant effects can generally be neglected if the 
Rabi frequency 



^2™ — 2Ai n d>n 



(8) 



is considerably less than the detuning of the field from off- 
resonant transitions. Although the situation is somewhat 
more complicated for fields with varying amplitude, off- 
resonant effects should remain small as long as 



max [2A m (t)d 



< Mm - Mn, Vn ^ TO, 



(9) 



i.e., as long as the peak amplitude 2A m of the control field 
f m is sufficiently small. Since the detuning of adjacent 
transitions for a Morse oscillator with anharmonicity a is 
fi m — fi m +i — 2a, the peak amplitudes of control pulses 
should be chosen such that 



A m <^ ct/d 

n 



(10) 
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The evolution of the controlled system is determined 
by a time-evolution operator U(t), which satisfies the 
Schrodinger equation 



i-U(t) = \H +J2H m [f m (t)]\u(t), (11) 



as well as the initial condition U(0) = I, where I is the 
identity operator. Note that throughout this paper we 
shall choose units such that U = 1. Our aim is to achieve 
a certain evolution of the system by applying a sequence 
of simple control pulses. Concretely, we seek to dynam- 
ically realize a desired unitary operator U at a certain 
target time T, i.e., 



U = U(T). 



(12) 



In some cases, we may not wish to specify a target time 
in advance, in which case we simply attempt to achieve 
the control objective at some later time T > 0. 



III. DETERMINATION OF PULSE SEQUENCE 
USING LIE GROUP DECOMPOSITIONS 

To solve the problem of finding the right sequence of 
control pulses, we apply the interaction picture decom- 
position of the time-evolution operator U(t), 



U(t) = U (t)Ui(t), 



(13) 



where Uo(t) is the time-evolution operator of the unper- 
turbed system, 



JY 



U (t) = cxp (-Mot) = J2 e ~ iEnt \ n )(n\ ( 14 ) 

71=1 

and Ui(t) comprises the interaction with the control 
fields. Inserting (|l^) into the Schrodinger equation (pT 
gives 

M 



H U(t) + H m [f m (t)]Uo(t)Ui(t), 



which leads to 



M 



ijUi{t) = U (rf \ H m [f m (t)] \ U (t)Ui(t). (15) 



Inserting equations ([l4| ) and (Q) into the right hand side 
of equation ( [l5"| ) leads after some simplification (see ap- 
pendix |a|) to 
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where we define 
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Hence, if we apply a control pulse 

fk(t) = 2A k (t) cos(n m t + <j) k ) 

= A k {t) +e -i(M m t+**)l ) (i 8 ) 

which is resonant with the transition frequency /x m , for 
a time period t k -i < t < tk and no other control fields 
are applied during this time period then we have 



U I {t) = V k (t)Ui(tk- 1 



(19) 



where the operator V k {t) is 



V k (t) = cxp 



A k (r) dr d m (sm<f> k 



>4>ky m ) 



(20) 

Thus, if we partition the time interval [0, T] into K subin- 
tervals [t k -\, t k ) such that to — and tx — T and apply a 
sequence of fixed-frequency control pulses, each resonant 
with one transition frequency /i m = ^J, a ( k ) , such that dur- 
ing each time interval only one control field is applied, 
then we have 

U(T) = U (T)Ui(T) = e- lkaT V K V K -i ■■■V 1 , (21) 
where the factors V k are 



% = exp 



A k (t) dr d„(u) (sin 4>kX<j{k) 



COS( 



'kya(k)) 



(22) 

2A k (t) is the envelope of the fcth pulse and it is a mapping 
from the index set { 1 , . . . , K} to the control index set 
{1, . . . , M} that determines which of the control fields is 
active for t € [tjt— ljifc]. 

It has been shown in |l5[| that any unitary operator 

U can be decomposed into a product of operators of the 
type 14 and a phase factor e* r = det U, i.e., there exists 
a positive real number T and real numbers C k and 4> k 
for 1 < k < K and a mapping a from the index set 
{1, . . . , K} to the control-sources index set {1, . . . , M} 
such that 



U = e* r V K V K _ 1 ...V 1 , 
where the factors are 



Vfc = exp [Cfc(sin< 



)] 



(23) 



(24) 



This decomposition of the target operator into a prod- 
uct of generators of the dynamical Lie group determines 
the sequence in which the lasers are to be turned on and 
off. A general algorithm to determine the Lie group de- 
composition for an arbitrary operator U is described in 
appendix M. 
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Note that in many cases the target operator U is 
unique only up to phase factors, i.e., two unitary op- 
erators U\ and U2 in U(N) are equivalent if there exist 
values 9 n € [0, 2ir] for 1 < n < N such that 



/ n n 
U2 = Ui £>""|n)<n| 

\n=l / 



(25) 



where |n) are the energy eigenstates. For instance, if 
the initial state of the system is an arbitrary ensemble of 
energy eigenstates 



N 

Po = ^2w n \n)(n\, 

n=l 



(26) 



where w n is the initial population of state \n) satisfying 



< W n < 1 and w n = 1, then we have 

JV 

U2poU% = Ux^\n)e iK w n e- iK {n\Ul 

n=l 



N 

= Ui^Wn\n){n\vt = UipoUi, (27) 

n=l 



i.e., the phase factors e 1011 cancel out. Thus, if the initial 
state of the system is an ensemble of energy eigenstates, 
which of course includes trivial ensembles such as pure 
energy eigenstates, then we only need to find a Lie group 
decomposition of the target operator U modulo phase 
factors, i.e., it suffices to find matrices V% such that 

u(yy^\n){n\\=V K V K _ l ...V 1 . (28) 



Note that decomposition modulo phase factors, when suf- 
ficient, is more efficient since it requires in general up to 
2(N — 1) fewer steps than the general decomposition al- 
gorithm. See appendix || for details. 



IV. AMPLITUDE AND PULSE LENGTH 

Comparing equations ( p2| ) and d^) shows that the con- 
stants C k in the decomposition (|25|J determine the pulse 
area of the fcth pulse: 



2A k (r)dT = 2 / A k {r)dT 



2C k 
d a {k) 



(29) 



There is considerable flexibility as regards the shape of 
the control pulses that are used. In this paper we consider 
two simple types of control functions: square- wave pulses 
with finite rise and decay times and Gaussian wavepack- 
ets. The former type of control pulses is convenient since 
it is easy to realize them the laboratory. In the opti- 
cal regime, for instance, a combination of CW lasers and 
Pockel cells can be used to achieve such pulse shapes. 



Gaussian control pulses can be derived from pulse laser 
systems. In addition, the latter type of control pulses 
has the distinct advantage of minimal frequency disper- 
sion, which is generally desirable for frequency-selective 
control. 



A. Square- wave pulses 

If we apply a field of constant amplitude 2Ak for a fixed 
period of time At k = t k — t k -\ then A k is determined by 
the pulse area constraint (|29|) 



A k = 



Ck 



(30) 



Thus, Ak can be adjusted by changing the pulse length 
Atfc, which allows us to account for laboratory con- 
straints on the strengths of the control fields and limit 
undesirable off-resonant effects by ensuring that 



A k < min \(i n -i - fi n \/d a 



(ft)- 



(31) 



Ideal square-wave pulses are rather simple and conve- 
nient pulse shapes. However, in practice it is impossible 
to reproduce square-wave pulses exactly. Rather, each 
pulse will have a certain rise and decay time To, which 
leads to pulse envelopes of the form depicted in Fig. 0. 
Mathematically, such pulses can be modeled by choosing 
the envelope of the pulse of the form 

2A k (t) = A fe {l + erf[4(<-ro/2)/r ]} 

+A k {1 + erf [4(i - At + r /2)/r ]} (32) 

where erf(x) is the error function 

2 



erf(x) 



dt. 



(33) 



Although this envelope function may appear quite com- 
plicated, it turns out that its pulse area is rather simple 
to compute. Observing that the area 
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FIG. 2: Square wave pulse of total pulse length Atk with rise 
and decay time to and peak amplitude 2A k 



since — erf(x) = erf(— a;), we note that the total pulse 
area of a modified square wave pulse of peak amplitude 
2Ak and total pulse length Atk > 2to with rise and decay 
time To is simply equal to the area of a rectangle of width 
Atk — To and height 2 Ak . Thus, the pulse area constraint 
( p9| ) leads to 2Ak{Atk — t ) = 2Ck/d a ( k ) or equivalently 



A,= 



C k 



rf CT (fe)(Atfe - T ) ' 



(34) 



Although this formula is very similar to (J30|), note the 
importance of including the rise and decay time of the 
pulse. Neglecting To amounts to overestimating the pulse 
area, which will result in poor control. 

Unless restrictions on the pulse strength dictate oth- 
erwise, the pulse lengths Atk can be chosen to be the 
same. For instance, if the target time for achieving the 
control objective is T and the Lie group decomposition 
shows that K pulses are required to achieve the control 
objective, then we would usually set 



T 

At k = — , 1 < k < K. 



(35) 



However, recall that it is important to assure that equa- 
tion (^|) is satisfied, i.e., that 2Ck/(Atk — t ) is much 
smaller than the detuning of the pulse frequency from off- 
resonant transitions. For instance, for a Morse oscillator 
with anharmonicity a, the detuning of each pulse from 
off-resonant transitions is at least 2a and thus Atk — To 
should be much larger than Ck /a in this case. 



B. Gaussian wavepackets 

If we wish to use Gaussian wavepackets instead of 
square waves then we choose the envelope of the fcth con- 
trol pulse to be of the form 



2A k (t) = 2A k e~ q ^ t - At ^ 2 ~ tk - 1 '> 2 



which corresponds to a Gaussian wavepacket centered at 
t* k = tk-i + Atk/2 of peak amplitude 2Ak- The peak 
amplitude 2Ak is determined by the pulse area constraint 
(P9|) . Concretely, we have 



-4 t 



qkC k 



(37) 



provided that the fcth time interval Atk is large enough 
to justify the assumption 



e -q 2 k (t-At k /2-t k - 
'tfc-1 

Noting that 



\2 

l) dt 



+ 00 



2 2 

■« T dr 



qk ' 



At k /2 

e 

-At fc /2 



2.2 

dt 



Qk 



-erf(<7 fc At fe /2) 



and erf(2) = 0.995322 we see, for instance, that choosing 

(38) 



qk 



4 
~AT k 



guarantees that more than 99% of the pulse area of the 
fcth pulse is contained in the control interval [tfc— i,tfe]- 
Choosing qk = 6/ Atk would ensure that more than 
99.99% of the pulse area of the fcth pulse is contained 
in the kth control interval. Again, unless restrictions on 
the pulse strength dictate otherwise, we set 



AU 



T 



l< k < K, 



(39) 



where T is the target time and K is the number of pulses 
required. Note that the peak amplitude 2Ak of the fcth 
field is proportional to qk , which is in turn inversely pro- 
portional to Atk- Hence, we can again limit the peak 
amplitudes and thus off-resonant effects by choosing the 
pulse lengths Atk sufficiently large. In the following sec- 
tions we apply these results to various control problems. 



V. POPULATION TRANSFER 1) -> \N) FOR A 
N-LEVEL SYSTEM 

Let us first consider a A^-level system initially in the 
ground state |1) and apply the decomposition technique 
to the simple problem of transferring the population of 
state |1) to the highest excited state \N) by applying a 
sequence of monochromatic control pulses, each resonant 
with one of the transitions frequencies n m , which can 
be regarded as the population inversion route to short- 
wavelength lasers. It can easily be verified that any evo- 
lution operator U of the form 



U 






A N -i ' 








(40) 



where A^-i is an arbitrary unitary (N — 1) x (N — 1) 
(36) matrix, e l9 is an arbitrary phase factor and is a N — 1 
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2A, 




FIG. 3: Gaussian wavepacket with qt = with peak am- 
plitude 2A k 



vector whose elements are 0, achieves the control objec- 
tive since 






A N -i 


e ie 








and thus the population of state \N) is equal to 
V e~ i0N e t0N = 1 after application of U. Next we observe 
that setting 

U = U {T)U I , Ui = V N -xV N -i->-Vx, (41) 
where the factors are 



V m 



exp 



- (sin0 m x m - coscj) m y m ) 



(42) 



) + £ 



for 1 < m < N — 1, always leads to a U of the form 
(pp|), independent of the initial pulse phases (f> m . This 
factorization corresponds to a sequence of N — 1 control 
pulses, where the mth pulse is resonant with the tran- 
sition | to) — > | to + 1) and has pulse area 7r/c? m , which 
transfers the population step by step to the target level 



N. 



field plot indicate that the corresponding pulse is reso- 
nant with the transition \m) — > \m + 1). The middle 
graphs shows the evolution of the energy-level popula- 
tions resulting from the application of the control fields, 
and the bottom graphs the energy of the system as a 
function of time, where the upper and lower horizontal 
lines in the bottom graph indicate the kinematical upper 
and lower bounds for the energy. Observe that for both 
choices of the control pulses the final population of tar- 
get level four is one, i.e., complete population transfer is 
achieved, while the populations of the intermediate levels 
increase and decay intermittently. 

Also note that the energy of the system increases 
monotonically from its kinematical minimum value at 
t = to its maximum value at the final time in both 
cases. The gradient of approach, however, is more uni- 
form for square-wave pulses, while Gaussian pulses tend 
to result in short, steep increases with long intermittent 
plateau regions. This is an advantage of using square- 
wave pulses, which could be significant especially for re- 
alistic systems where spontaneous emission is a prob- 
lem. Gaussian wavepackets on the other hand, have the 
distinct advantage of minimal frequency dispersion and 
should thus be less likely to induce unwanted off-resonant 
effects. 

As regards off-resonant effects in general, note that 
Cfe = 7r/2 for all pulses. Thus, for our Morse oscillator 
with a — 0.1 and square- wave pulses with rise and decay 
time To = 30 time units, we have Ck/ot = 5tt, which is 
much smaller than our choice of Atk — To = 170 time 
units. More precisely, the peak amplitude of field f m , 
which is resonant with the transition |m) — > \m + 1), is 



2A r . 



7T TiLOo 

(200- 30)^ PiT' 



Recalling that the dipole moment is d m = \pm pi2, the 
"Rabi-frequency" of the pulse is thus 2A m d m /T% = 7r/170 
wo 7 which is less than 1/10 of 0.2 ui , the detuning of 
the pulse frequency from closest off- resonant transitions. 
Note that for fixed pulse lengths, the peak amplitudes 
for Gaussian control pulses are necessarily larger than for 
square- wave pulses, in our case by about a factor of 1.5. 
However, it seems reasonable to assume that Gaussian 
control pulses that have the same pulse area and pulse 
length as square- wave pulses should be less likely to cause 
strong off-resonant effects than square- wave pulses due to 
smaller frequency dispersion. 



Thus, the control sequence derived using the Lie group 
decomposition technique agrees in this case with the ob- 
vious choice of the control pulses. 

Figure |J shows the results of control computations 
for the four-level Morse oscillator described above using 
square- wave and Gaussian control pulses, respectively. 
The top graph in each figure shows the envelopes of the 
control pulses in units of 10 _2 ?iwo/pi2j where w is the 
oscillator frequency and p 12 is the dipole moment of the 
1 — ► 2 transition. The labels f m for m = 1,2, 3 in the 



VI. INVERSION OF ENSEMBLE 
POPULATIONS FOR A MIXED-STATE SYSTEM 

The results of the previous section are encouraging 
in that they agree with intuitive control schemes. The 
power of the decomposition technique, however, lies in 
its ability to predict control schemes for problems that 
have no such obvious solutions. One such case is a gen- 
eralization the control problem discussed in the previous 
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FIG. 4: Population transfer from the ground state |1) to the excited state |4) for the four-level Morse oscillator described above 
using three square-wave control pulses with rise and decay time to = 30 time units (a) and Gaussian pulses with shape factor 
q = 4/100 (b). The mth pulse is resonant with the transition \m) —* \m + 1} and has pulse area n/d m for m = 1, 2, 3. 



section, where we attempt to achieve complete inversion 
of the ensemble populations for a TV-level system whose 
initial state po is an arbitrary ensemble of energy eigen- 
states (p^). This control operation can be regarded as 
a form of generalized NOT gate for mixed-state iV-level 
systems, the ensemble NOT gate, which should not be 
confused with the U-NOT gate described by Buzek et. al. 
in HI. 



reversal of the ensemble populations for the system is 



U = 



( • • • e Wl \ 



\ 







e l 








0/ 



(43) 



The desired evolution operator to achieve a complete 



where the e l9n are arbitrary phase factors. Assuming 
as above that each transition between adjacent energy 
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levels can be individually addressed by selecting the fre- 
quency of the control pulse and possibly by using other 
selection rules, the generators of the dynamical Lie alge- 
bra are again of the form (|22] ) and a possible Lie group 
decomposition in terms of these generators is given by 



u=mt) n 



e=N-i 



(44) 



Also note that the energy of the system increases 
monotonically from its minimum value in thermal equi- 
librium to its kinematical maximum value at the final 
time, as expected, where the gradient of approach is 
more uniform for square-wave pulses than for Gaussian 
wavepackets. As regards off-resonant effects, the same 
considerations as in the previous section apply, since 
Cfc = 7r/2 for all pulses, as in the previous example. 



where the factors V m are as defined in ( |42|) . This decom- 
position corresponds to a sequence of K = N(N — l)/2 
control pulses 



h,f2 



, /iV-1, fl,h, ■ ■ ■ , /iV-2, . . . ,/l,/2, fl, 



where the fcth pulse is resonant with the transition 
|cr(fc)) — > \cr(k) + 1) and has pulse area Tv/d^^y This 
decomposition is optimal in the sense that a complete in- 
version of the ensemble populations cannot be achieved 
with K' < K control pulses if the initial populations w n 
satisfy w n ^ w m for n 7^ m. 

As an example, we again apply the scheme to the N- 
level Morse oscillator with energy levels (^|) and transition 
dipole moments (||) described above. To be specific, we 
assume that the system is initially in thermal equilibrium 
at a temperature (E^ — Eij/k, i.e., the initial populations 
of the energy eigenstates are given by the Boltzmann 
distribution 



Wn 



exp[(E n - E^/jEN - E,)} 
J2k^exp[(E k - E 1 )/(E N - Et)} 



(45) 



for n — 1,2,3,4. Note that the initial populations sat- 
isfy w\ < W2 < W3 < W4. Our goal is to create an 
anti-thermal ensemble, i.e., an ensemble where the popu- 
lations of the energy eigenstates are reversed so that the 
ground state |1) has the lowest population W4, and the 
highest excited state has the highest population wi, etc. 
Figure [| shows the result of control simulations using 
square-wave and Gaussian controls, respectively. Note 
that each pulse in the control sequence interchanges the 
populations of two adjacent energy levels until a complete 
reversal of the populations is achieved; for our four-level 
system with initial populations w n the effect of the con- 
trols on the populations can be summarized as follows 







fl 


h 




h fl 




h 


fl 


|1>: 


Wi 


> 


W 2 — > 


w 2 


-> W 2 -* 


w 3 


-> w 3 


—y W4 


|2>: 


W 2 




Wi -> 


U> 3 


-> W 3 


W2 


— > W4 


-> w 3 


|3>: 


W 3 




W 3 -> 


Wi 


— > W4 — * 


W4 


— > w 2 


— > W2 


|4): 


W4 




W4 — » 




-> Wi -> 




-> Wl 





where / m refers to a control pulse of frequency jj, m and 
pulse area 7r/rf m for m = 1,2,3. The first pulse inter- 
changes the populations of levels |1) and |2), the second 
pulse flips the populations of levels |2) and |3), the third 
pulse switches the populations of levels |3) and |4), etc. 



VII. CREATION OF ARBITRARY 
SUPERPOSITION STATES 

In this section we consider the problem of creating ar- 
bitrary superposition states from an initial state. Con- 
trol schemes that create such superposition states may 
prove very useful in controlling quantum interference in 
multi-state systems and can be considered a generaliza- 
tion of 7r/2 pulses used routinely in free induction-decay 
experiments. To be specific, we assume that the sys- 
tem is initially in the ground state |1) and that the goal 
is to achieve a superposition state \^(t)). Observe that 
any (normalized) superposition state can be conveniently 
written as 



!*(*)) 



N 

E 

71=1 



r n e 



»(<))» 



(46) 



where \n(t)) — e lEnt \n) is a rotating frame and the coef- 
ficients r n satisfy the normalization condition Xm=i r n 
1 . Thus, in order to reach the target state at time 

T we need to find a unitary operator U(T) such that 



U(T) 







r±e 
r 2 e 



i(8i 
i(S 2 



-E t T) 
-E 2 T) 



\ 



\ r N e 



i(9 N -E N T) 



(47) 



However, note that it actually suffices to find a unitary 
operator U\ such that 



Ui 







( n \ 

''2 



(48) 



V°/ VnJ 

since if Ui satisfies © then U(T) = ®{T)Ui where 

9(T) = dmg{e l ^- E ' T \ . . .^-e^ (4Q) 

automatically satisfies (|4^). In order to find a unitary 
operator U\ that satisfies dig) we set 

/ n \ 

C/(°)= r " . (50) 



\ r N 



lN-1 



J 
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FIG. 5: Population inversion for a four-level Morse oscillator initially in thermal equilibrium using six square-wave control 
pulses with rise and decay time to = 30 time units (a) and Gaussian pulses with shape factor q — 4/100 (b). Each pulse labeled 
fm is resonant with the transition \m) — > \m + 1) and has pulse area n/d m . 



where jjv-i is the identity matrix in dimension N — 1, the superposition state 
and perform Gram-Schmidt orthonormalization on the 
columns of . The resulting matrix U\ will be uni- 
tary and satisfy JIs| ) . To determine the control pulse se- 
quence, we then apply the decomposition algorithm pre- 4 ^ 4 
sented in appendix [B]to factor the target operator \^(t)) — - \n(t)) = - ^ e~ iEnt \n) (52) 



2 

n=l 



U = Uo(T)^(T)U 1 = Q(Qpi. (51) 



As an example, let us consider the problem of creating 



for a four-level system initially in state |1). In order to 
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find a unitary operator U such that 

(A (w\ 

1/2 
1/2 



U 



(53) 



V 1 / 2 / 



we set 



( 1/2 0\ 

1/2 10 

1/2 10 

V 1/2 o o i y 



(54) 



and perform Gram-Schmidt orthonormalization on the 
columns of U^°\ which gives 



( 1/2 -V3/6 -V6/6 -V2/2\ 

1/2 +V3/2 

1/2 -V3/6 +V6/3 

V 1/2 -V3/6 -a/6/6 +a/2/2 / 



(55) 



Since O(T) = Uq(T), the target operator is simply 

the decomposition al- 
leads to the following 



U (Tye(T)U = U x . Applying 
gorithm described in appendix n 
factorization olU\. 



U x = VsVaVMVx, 
where the factors are 



(56) 



V\ = exp(+Cia;i) 

V 2 = exp(-C 2 x 2 ) 

V 3 = cxp(+C 3 i 3 ) 

V4 = CXp ( + C4X2) 

V5 = exp(-C 5 xi) 



1 - 3> 

C 2 = arctan (V%) 

7T 

3—4! 
(~i 7T 

O4 - 2' 
5—9- 



(57) 



This decomposition corresponds to the following se- 
quence of five control pulses 



hit) 


= A 1 (t)e i ^ t+7r ^ - 


- c.c. 


= -2Ax(t) sin(/iit) 


hit) 


= A 2 {t)e l ^ t - 7T IV - 


- c.c. 


= +2A 2 (t) sin(^ 2 i) 


hit) 


= A 3 {t)e l ^ t+7r ^ - 


- c.c. 


= -2A 3 (t) sin^ 3 t) 


hit) 


= A 4 (t)e l ^ 2t+7r /V _ 


1- c.c. 


= -2A A (t) sin( ( u 2 i) 


hit) 


= Ag^e 1 ^ 1 *-^ 2 ) - 


- c.c. 


= +2A 5 (t) sm(fiit) 



with pulse areas 2n/3di, 2 arctan(v / 2)/rf2, n/2d 3 , n/d 2 
and 7r/c?i, respectively. Note that only five instead of 
six pulses are required in this case since the target op- 
erator U\ has two consecutive zeros in the last column, 
which implies that one of the six control pulses has zero 
amplitude and can thus be omitted. The results of two 
control simulations based on this decomposition of Ux 
using square-wave and Gaussian control pulses, respec- 
tively, are presented in figure ^. Observe that the expec- 
tation value of the projector onto the target state 
assumes its upper bound of one at the final time in both 



cases, which implies that the system has indeed reached 
the target state 

Unlike decompositions ( |4l| ) and ([ll]) in the previous 
sections, in which the initial phases (f> m of the control 
pulses were arbitrary, the factorization ( [56] ) fixes the 
pulse area and frequency [i m as well as the initial phase 
(f> m of each pulse. While the dependence of the control 
scheme on the pulse area and frequency is expected, the 
dependence on the initial phases of pulses may be a rea- 
son for concern. Thus, the question arises how the initial 
phases of the control pulses affect the outcome of the 
control process. In order to answer this question, we 
compute the unitary operator 



U 2 = V 5 V 4 V 3 V 2 V 1 , 
where the factors are 



(58) 



Vi 


= exp[Ci(sin( 


hx\ 


— COS ( 


km)} , 




= exp [C 2 (sin ( 


fox 2 


— COS< 


fam)] > 


v 3 


= exp [C 3 (sin < 


hx 3 


— COS< 


fafc)] , 


v 4 


= cxp[C4(sin( 


hx 2 


— COS ( 


hfa)} , 


% 


= exp [C5 (sin < 


kxi 


— COS ( 


hm)] 



(59) 



and the constants Ck are as in (|57|) but the initial phases 
4>k of the control pulses are allowed to be arbitrary, and 
apply this operator to the initial state |1). The resulting 
state 






W 



/ \/2 e «(04 + </>5-01-02) \ 

1/2 e i(-T/2-<fe) 

1/2 eK 7 ^^ 1- ^ 4 ) 
\ 1/2 eK^/ 2 ^ 1- * 2- * 3 ) J 



(60) 



differs from the desired target state only in the phase 
factors, i.e., the pulse phases do not affect the relative 
amplitudes r n of the superposition state created. Fur- 
thermore, we can make the phase factors equal one by 
choosing 4> 2 = n/2 — 2<f>x, (f> 3 = <j>i, <f>4 = n — 4>\ and 
05 = — 7r/2, where <f>\ can be arbitrary. Note that the 
choice 4>i = 7r/2, (j) 2 = —n/2, <f> 3 = tt/2, $4 = tt/2 and 
05 = —n/2 we made above satisfies these conditions. Re- 
calling that we expanded the state with respect to a ro- 
tating frame \h(t)) — e~ lEnt \n), shows that the phases 
associated with the energy eigenstates in this superposi- 
tion state are oscillating rapidly compared to the length 
of the control pulses. It is therefore hardly surprising 
that control of relative phases between the energy eigen- 
states in the superposition state requires control over the 
pulse phases. 



VIII. OPTIMIZATION OF OBSERVABLES 

Finally, we address the problem of maximizing the en- 
semble average of an arbitrary observable for a system 
whose initial state is an ensemble of energy eigenstates 
(p6|) . Let us first consider the case of a time-independent 
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FIG. 6: Creation of the superposition state = \ y^ 4 _. |w(t)) for a four-level Morse oscillator initially in the ground state 

|1) using square- wave control pulses with rise and decay time to = 30 time units (a) and Gaussian control pulses with shape 
factor q = 4/100 (b). Each pulse labeled f m is resonant with the transition \m) — > \rn + 1). The pulse areas are 2n/3di, 
2 arctan(\/2)/rf2, Tv/2d,3, iv/d2 and n/di, respectively. 



observable A. To determine the target operator required 
to maximize the ensemble average (A) of A we observe 
that (A) is bounded above by the kinematical upper 
bound Il8| 



where A„ are the eigenvalues of A counted with multi- 
plicity and ordered in a non-increasing sequence 



Ai > A 2 > • • • > A 



JV- 



(62) 



JV 



n=l 



(61) 



and w n are the populations of the energy levels E n of 
the initial ensemble and a is a permutation of {1, ... , N} 
such that 



UV(1) > W a {2) >■> Wa(N)- 



(63) 
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Observe that this universal upper bound for the ensemble 
average of any observable A is dynamically attainable 
since the systems considered in this paper are completely 
controllable |]| ||. 

Let l^n) for 1 < n < N denote the normalized eigen- 
states of A satisfying A|^„) = \ n \^> n ) and let Ux be a 
unitary transformation such that 

\* a(n) ) = Ui\n), l<n<N. (64) 

Given an initial state po of the form (p6|), we have 



Tr (AUxpoUl 



Tr \ AJ2wnUi\n){n\Ulj 
Tr ^u; n i|* CT(n) )(* ff(n) |^ 

Tr w n^a(n)\^ a(n)) crin)^ 

^W n A CT ( n ) = ^ W <r(n)^n- (65) 



Hence, if the system is initially in the state ( p!6"| ) then XJ\ is 
a target operator for which the observable A assumes its 
kinematical maximum and we can use the decomposition 
algorithm described in appendix |b| to obtain the required 
factorization of the operator Uj = Uq(T)<Ui. 

However, if A is an observable whose eigenstates are 
not energy eigenstates, then the expectation value or en- 
semble average of A will usually oscillate rapidly as a 
result of the action of the free evolution operator Uo(t). 
These oscillations are rarely significant for the applica- 
tion at hand and often distracting. In such cases it is 
advantageous to define a dynamic observable 



A(t) = U (t)AU (tf 



(66) 



and optimize its ensemble average instead. To accom- 
plish this, note that if |\I> n ) are the eigenstates of A sat- 
isfying i|*„) = \ n \$ n ) then |#n0)) = U(t)\^ n ) are the 
corresponding eigenstates of A(t) since 

A(t)\V n (t)) = U{t)AU{tfU{t)\y n ) 
= U(t)\ n \V n ) 
= A„|# n (i)) 

and if Ux is a unitary operator such that equation (|64| ) 
holds then Uo(t)Ui is a unitary operator that maps the 
energy eigenstates |n) onto the A(t)-eigenstates |^ n (t)) 
since 

U (t)Ui\n) = 0b(t) !¥«,(»)) - I 

for 1 < n < N. Thus, the evolution operator required 
to maximize the ensemble average of A(t) at time T > 



is Uo(T)Ui and the target operator to be decomposed as 
described in appendix O is 



U = U Q (T^U (T)Ux = Ux. 



(67) 



For instance, suppose we wish to maximize the en- 
semble average of the transition dipole moment operator 
A(t) = U (t)AU Q (tf, where 



N-l 



A = ^d n {\n){n + l\ + \n + l){n\) 



for a system initially in state ( |26| ) with 
Wi > W 2 > ■ ■ ■ > wn > 0. 



(68) 



(69) 



Then we need to find a unitary operator Ui that maps 
the initial state |n) onto the A-eigenstate \^ n ) for 1 < 
n < N. To accomplish this we let Ux be the N x N 
matrix whose nth column is the normalized A-eigenstate 
|* n ). Then U\ clearly satisfies Ui\n) = |^„) and since 
| * m ) = § ntm by hypothesis, i.e., the eigenstates are 
orthonormal, U\ is automatically unitary. 

For N = 4 and d n — y/n the eigenvalues of the operator 
A defined in equation (|6^) are (in decreasing order) 



Ai = y3 + V6, A 2 = V 3 - V6, A 3 = -A 2 , A 4 = -Ax 

and the corresponding eigenstates with respect to the 
standard basis \n) are the columns of the operator 



Ux 



1 


1 


1 


1 


2Ai 


2A 2 


2A 2 


2Ai 


1 


1 


1 


1 


2 


2 


2 


2 


V2+V3 


V2-V3 


V2-V3 


V2+V3 


2Ai 


2X 2 


2A 2 


2Ai 


1 


1 


1 


1 


2 


2 


2 


2 



(70) 



Applying the decomposition algorithm described in ap- 
pendix p| yields the following factorization 



where the factors are 


fi 


= exp (- 


-Cxxx), 


% 


= exp (- 


-C 2 x 2 ) , 


% 


= exp (- 


-C 3 xx), 


v A 


= exp (- 


-C4X3) , 


% 


= exp (- 


-C 5 x 2 ) , 


% 


= exp (- 


-C 6 xx) , 



Ux& = V 6 V 5 V A V 3 V 2 Vx 



Cx = tt/4, 

C 2 = arctan (V2) , 

C 3 = arccot ^-v|+3^ 

d = vr/3, 



(71) 



(72) 



C 5 = arctan 

C5 = arccot fV3 + \/6 



and 6 = diag(l, -1,1,-1). Note that U 2 = UxQ is 
equivalent to U\ since O commutes with po as defined 
in equation (|2q), i.e., Qpo& = po, and thus 



Tr 



(Au 2 p U^ = Tr (AUxQpo&Ul) 

= Tr (AUxPoUx) ■ (73) 
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Decomposition ( |71| ) corres 
qucnce of six control pulses 



sponds to the following se- 



hit) 


= Aiity^ 1 *- 


-T/2) - 


- C.C. 


= 2A 1 (i) sm(fiit) 


hit) 


= A 2 it)e l( ^ t ' 


-T/2) _ 


- C.C. 


= 2A 2 (t) sm(fi 2 t) 


hit) 


= A 3 (i)e i(Mlt " 


-w/2) . 


- C.C. 


= 2A 3 (t) sin(Mit) 


hit) 


= A^ty^ 1 - 


-T/2) . 


- C.C. 


= 2A 4 (t) sin(M3*) 


hit) 


= A 5 (t)e l ^ 2t - 


-w/2) . 


- C.C. 


= 2A 5 (i) sin(// 2 t) 


hit) 


= Aeity^ 1 *- 


-t/2) - 


- C.C. 


= 2A 6 (t) sin(^ii) 



with pulse areas ■n/2d\, 2C 2 jd 2 , 2C 3 /d\, 2-K/2>dz, 2G^jdi-, 
2Ce/di, respectively. Again, the decomposition fixes the 
frequency and pulse area as well as the initial phase of 
each pulse, and the question thus arises, what role the 
phases play. As we have already seen above, the target 
operator Ui is not unique. In fact, equation (|7^) shows 
that right multiplication of U\ by any unitary matrix 
that commutes with po produces another unitary oper- 
ator that is equivalent to Ui in that both evolution op- 
erators lead to the same ensemble average of the target 
observable. Nevertheless, in general, the control process 
is sensitive to the phases 4> m . For instance, one can verify 
that changing the phase <\>\ of the first pulse from — 7r/2 
to 7r/2 in the pulse sequence above leads to the following 
evolution operator 



U 3 
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2A 2 


2A X 


2A 2 


2A X 


1 
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V2-V3 


V2+V3 


V2-V3 




2X 2 


2Ai 


2A 2 


-2Ai 
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2 


2 


2 



(74) 



which maps |3) onto 1^3) and |4) onto — ^4} but |1) onto 
\ i &2) and |2) onto l^i) and thus leads to the ensemble 
average 



(A) = w\\ 2 + W2A1 + W3A 3 + W4A4 



(75) 



at the final time, which is strictly less than the kinemat- 
ical maximimum if w\ > u>2- 

Figure fl] shows the results of control simulations based 
on the control pulse sequence above for our four-level 
Morse oscillator initially in thermal equilibrium using 
square- wave and Gaussian control pulses, respectively. 
Observe that the observable indeed attains its kinemati- 
cal upper bound at the final time, as desired. 



IX. CONCLUSION 



We have presented several control schemes designed to 
achieve a wide variety of control objectives ranging from 
population transfers and inversions of ensemble popula- 
tions to creation of superposition states and optimization 
of observables. The main advantage of the schemes is 
their simplicity: the control objective is achieved by ap- 
plying a sequence of simple control pulses such as square- 
wave pulses with finite rise and decay times, which can be 
created in the laboratory by modulating the amplitude of 
CW lasers using Pockel cells, or Gaussian wavepackets, 
which can be obtained from pulsed laser sources. The 
main shortcoming of this control approach is the need to 
be able to address a single transition at a time, which 
limits the applicability of this technique to systems with 
sufficiently distinct transition frequencies (as regards the 
controlled transitions) to permit frequency-selective con- 
trol, unless other selection rules can be employed to ad- 
dress transitions individually. Another possible compli- 
cation arises from unwanted off-resonant effects. How- 
ever, as we have demonstrated, in general it is possible 
to limit such effects by using sufficiently weak control 
pulses. 



APPENDIX A: DERIVATION OF EQUATION @ 

Inserting equations (|lj) and (0) into equation ( |l5| ) and recalling [i m = E m+ i — E m leads to 



.d&iit) 
dt 
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EJE n t~ (a (f\J{V>mt+$m) J p ,1 4- A (t V~* Cm ™ t +<f>m > fl P ^ p-* S n'*S , ,fj T (t\ 

c c n,n 1 V W c ^m^m^rn-]-! \ -^m \ ° J Uj m c 'm-\-l,m J c c n ,n u 1 V ) 
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FIG. 7: Maximization of the (dynamic) transition dipole moment operator A(t) for a four-level Morse oscillator initially in 
thermal equilibrium using square-wave control pulses with rise and decay time to = 30 time units (a) and Gaussian control 
pulses with shape factor q = 4/100 (b). 



)Uj(t). 



Hence, multiplying both sides by — i gives 

dUr(t) _ 
dt ~ 



in 

(sm<p m x m — cos (j) m y m ) U i (t). 



(Al) 
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APPENDIX B: LIE GROUP DECOMPOSITION ALGORITHM 



To find a decomposition of the form (|23|) for a given 
unitary operator U, we set 



£(o) = e -;r/iV£ 



(Bl) 



where e ir = det(U), to ensure that S SU{N). Our 
goal is to reduce this matrix EA°) step by step to a diag- 
onal matrix whose diagonal elements are arbitrary phase 
factors e l6n . Recall that this reduction is always suffi- 
cient if the initial state of the system is an ensemble of 
energy eigenstates. 

Let lljj denote the zth row and jth column entry in 

the matrix representation of . In the first step of the 
decomposition we wish to find a matrix 



= exp [— Ci (sin — cos^iyi)] , 



(B2) 



which is the identity matrix everywhere except for a 2 x 2 
block of the form 



( cos(C*i) ie^ 1 sin(Ci) 
y ie'^ 1 sin(Ci) cos(C*i) 

in the top left corner, such that 



WW 



where c is some complex number. Noting that 



(B3) 



rj(0) 






V '■ J 







(B4) 



MO) 

U 1,N 



rie 



U. 



(0) 



2,N 



r 2 e 



it can easily be verified that setting 



<f>k = tt/2 + a% - a 2 
C k = -arccot(-r 2 /ri) 



achieves (B4). Next we set 

and find of the form 

VF (2) = exp [~C 2 (sin (f> 2 x 2 - cos 4> 2 y 2 )} 
such that 



(B5) 

(B6) 
(B7) 



/ \ 



u. 



(1) 



u. 



2, N 

(1) 

3, N 



V 



/ 





c 



(B8) 



where c is again some complex number. Repeating this 
procedure N — 1 times leads to a matrix ^( JV_1 ) whose 









exp 




x exp 



last column is (0, . . . , 0, e ieN ) T . Since we are not con- 
cerned about the phase factor e N in this paper, we stop 
here. Noting that 

r 7T 

- (sin0xjv_i - cos#at_i) 

with 4> = -7r/2 - 6<„ maps (O^* 9 "- 1 ) 11 onto (0, 1) T , we 
see that complete reduction to the identity matrix would 
require two additional steps at this point to eliminate the 
phase factor e l6,JV , which would result in two additional 
control pulses. 

Having reduced the last column, we continue with the 
(N — l)st column in the same fashion, noting that (since 
is unitary) at most N — 2 steps will be required 
to reduce the (N - l)st column to (0, . . . , 0, e i0N -\ 0) T . 
We repeat this procedure until after (at most) K = 
N(N - l)/2 steps C/ (0) is reduced to a diagonal matrix 



diag(e lSl , . 



r ), as required and we have 



,,e 



(B9) 



.. .W/W^o) =d iag(e l9 \ 
Finally, setting 

V k = (#< A '+ 1 - fc )) t (B10) 

leads to 

U {0) = VkVk-i ■ ■ ■ Vidiag {e l8 \ . . . ,e w ») (Bll) 
and therefore 

u = v K v K _ 1 ---v 1 e, 

where Q is a diagonal matrix 
9 = e^diag (e 1 



(B12) 
(B13) 



Note that U can always be decomposed such that Q is 
the identity matrix. However, to achieve this goal up to 
2(N — 1) additional terms would be required to eliminate 
the phase factors, which would result in additional 
control pulses. While some applications may indeed 
require the elimination of these phase factors, the phase 
factors are often insignificant and the additional control 
pulses would be superfluous. For a more sophisticated 
decomposition algorithm that requires only very few 
phases the reader is referred to |2l| . 
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